Resources:
There are 2 ways to install sgpv package.
install.packages("sgpv")devtools::install_github("weltybiostat/sgpv")#GitHub
install.packages("devtools")
devtools::install_github("weltybiostat/sgpv")
#CRAN
install.packages("sgpv")
?sgpv::sgpvalue
Functions Included:
sgpvalue()
sgpower()
plotsgpv()
plotman()
fdrisk()
#Setup
theta0 = 0
point = 0.49
data.lo = 0.15
data.hi = 0.82
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -0.125,
null.hi = 0.125)
## $p.delta
## [1] 0
##
## $delta.gap
## [1] 0.2
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -0.2083,
null.hi = 0.2083)
## $p.delta
## [1] 0.08701493
##
## $delta.gap
## [1] NA
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -0.25,
null.hi = 0.25)
## $p.delta
## [1] 0.1492537
##
## $delta.gap
## [1] NA
sgpvalue(est.lo =data.lo,
est.hi = data.hi,
null.lo = -0.5,
null.hi = 0.5)
## $p.delta
## [1] 0.5223881
##
## $delta.gap
## [1] NA
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -0.75,
null.hi = 0.75)
## $p.delta
## [1] 0.8955224
##
## $delta.gap
## [1] NA
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -1,
null.hi = 1)
## $p.delta
## [1] 1
##
## $delta.gap
## [1] NA
#SBP Example Background
# Data means
xbar=c(141,142,143.5,144,146,145,145.5,146)-1
# Data Standard Errors
se=c(.5,1,.5,1,2.25,1.25,.25,.5)
# Indifference Zone
delta.a=143
delta.b=147
#Point Null
h0=145
# Estimated Uncertainty Interval Bounds
lb<-xbar-1.96*se
ub<-xbar+1.96*se
sgpvalue(est.lo=lb,
est.hi=ub,
null.lo = delta.a,
null.hi = delta.b)
## $p.delta
## [1] 0.0000000 0.0000000 0.2448980 0.5000000 0.5000000 0.7040816 1.0000000
## [8] 1.0000000
##
## $delta.gap
## [1] 1.01 0.02 NA NA NA NA NA NA
plotsgpv(est.lo = lb,
est.hi = ub,
null.lo = delta.a,
null.hi = delta.b,
plot.axis=c("TRUE","FALSE"),
null.pt=0, outline.zone=TRUE,
title.lab="SBP Example",
x.lab="Classical p-value ranking",
legend.on=FALSE)
# sgpower(true=xbar,
# null.lo=delta.a,
# null.hi=delta.b,
# std.err=se,
# interval.type='confidence',
# 'interval.level'=0.05)
# Raw p-value Z-value
z1<-(xbar-h0)/se
# Max p-value Z-value
z2<-ifelse(abs(xbar-delta.b)< abs(xbar-delta.a),(xbar-delta.b)/se,(xbar-delta.a)/se)
z2<-ifelse(xbar>delta.b | delta.a > xbar,z2,0)
# Raw p-value
p1 <- round(2*pnorm(-abs(z1)),4)
# Max p-value
p2 <- round(2*pnorm(-abs(z2)),4)
#SGPV
p3 <- sgpvalue(est.lo=lb,
est.hi=ub,
null.lo = delta.a,
null.hi = delta.b)$p.delta
## Plot of SBP example
id.color=rgb(208,216,232,max=255)
plot(c(xbar,NA),1:(length(xbar)+1),xlim=c(min(lb,delta.a),max(ub)+5),xaxt='n',yaxt='n',ylab="",xlab="",type="n")
rect(delta.a,length(xbar)+1,delta.b,0,col=id.color,border=NA)
abline(v=h0,lty=2,lwd=2,col='black')
points(c(xbar,NA),1:(length(xbar)+1), pch=20)
for(i in 1:8){
lines(c(lb[i],ub[i]),c(i,i))
}
text(max(ub)+4.5,1:(length(xbar)+1),c(round(p3,4),"SGPV"))
text(max(ub)+2.5,1:(length(xbar)+1),c(round(p2,4),"Max p-value"))
text(max(ub)+.5,1:(length(xbar)+1),c(round(p1,4),"p-value"))
axis(side = 1,at=seq(round(min(lb)),round(max(ub)),delta.b-h0))
axis(side=1,at=seq(144,148,2)) # to re-draw axis
mtext(side=1,"Systolic blood pressure (mmHg)",line=2.5)
###### BOMAMI Trial Slide 28 and 29
## OR
sgpvalue(est.lo=0.967,
est.hi=7.286,
null.lo=0.9,
null.hi=1.11)$p.delta
## [1] 0.3404762
## RR
sgpvalue(est.lo=0.953,
est.hi=4.589,
null.lo=0.9,
null.hi=1.11)$p.delta
## [1] 0.3738095
## RD
sgpvalue(est.lo=0.003,
est.hi=0.352,
null.lo=-0.1,
null.hi=0.1)$p.delta
## [1] 0.277937
## Logistic Regression Slide 30
###### BOMAMI Trial
## Primary OR
p1 = sgpvalue(est.lo=1.18,
est.hi=20.19,
null.lo=0.9,
null.hi=1.11)$p.delta
## Covariates
p2 = sgpvalue(est.lo=1.09,
est.hi=19.28,
null.lo=0.9,
null.hi=1.11)$p.delta
p3 = sgpvalue(est.lo=0.86,
est.hi=38.29,
null.lo=0.9,
null.hi=1.11)$p.delta
p4 = sgpvalue(est.lo=0.66,
est.hi=8.98,
null.lo=0.9,
null.hi=1.11)$p.delta
p5 = sgpvalue(est.lo=0.36,
est.hi=26.92,
null.lo=0.9,
null.hi=1.11)$p.delta
p6 = sgpvalue(est.lo=0.95,
est.hi=1.11,
null.lo=0.9,
null.hi=1.11)$p.delta
p7 = sgpvalue(est.lo=0.09,
est.hi=1.59,
null.lo=0.9,
null.hi=1.11)$p.delta
# Table with results
kable(cbind(c("Variable",
"group",
"tobacco",
"microvascular obstroction",
"dyslipidemia",
"gender",
"age",
"hypertension"),
c("Odds Ratio", 4.89,4.59,5.72,2.32,3.12,1.02,0.38),
c("Raw p-value", 0.03,0.04,0.07,0.22,0.30,0.56,0.19),
c("SGPV",p1,round(p2,3),p3,p4,p5,p6,p7)), format = "html", table.attr = "style='width:50%;'") %>%
kable_styling(c("striped", "bordered"))
| Variable | Odds Ratio | Raw p-value | SGPV |
| group | 4.89 | 0.03 | 0 |
| tobacco | 4.59 | 0.04 | 0.048 |
| microvascular obstroction | 5.72 | 0.07 | 0.5 |
| dyslipidemia | 2.32 | 0.22 | 0.5 |
| gender | 3.12 | 0.3 | 0.5 |
| age | 1.02 | 0.56 | 1 |
| hypertension | 0.38 | 0.19 | 0.5 |
Survival time in patients with advanced lung cancer (days)
Potential for gender dissimilarities
Trial by North Central Cancer Treatment Group (1994): https://pubmed.ncbi.nlm.nih.gov/8120560/
Interval Null [-0.05, 0.05] % difference
# Survival Model stratified females and Males
fit<- survfit(Surv(time, status) ~ sex, data = lung)
# plot(fit,col=c("red","blue"))
# summary(fit)
# objects(fit)
# summary(fit)$surv
newdata=lung$sex
lung.Surv <- with(lung, Surv(time=time, event=status))
lung.survfit <- survfit(lung.Surv ~ lung$sex)
sCox <- coxph(lung.Surv ~ as.factor(sex),data=lung)
#Compute SGPV of difference
pred.1=summary(survfit(sCox,newdata=data.frame(sex=1)))$surv
pred.2=summary(survfit(sCox,newdata=data.frame(sex=2)))$surv
pred.diff=pred.2-pred.1
time.diff=summary(survfit(sCox,newdata=data.frame(sex=1)))$time
v.1=summary(survfit(sCox,newdata=data.frame(sex=1)))$std.err^2
v.2=summary(survfit(sCox,newdata=data.frame(sex=2)))$std.err^2
se.diff=sqrt(v.1+v.2)
lb=pred.diff-1.96*se.diff
ub=pred.diff+1.96*se.diff
pdelt <- sgpvalue(lb,ub,-0.05,0.05)$p.delta
#Survival Plot
plot(lung.survfit,col=c("dodgerblue1","hotpink"),mark.time=F,lwd=c(2,2),
ylab="Survival",xlab="") ## sex=1 is hotpink
for (i in 1:length(time.diff)){rug(time.diff[i],col=COL[i],lwd=1.2,ticksize=0.04)}
axis(side=1)
mtext("Days",side=1,line=2.25)
legend(750,1,c("Females","Males"),col=c("hotpink","dodgerblue1"),lty=1,lwd=2,bty="n")
par(xpd=TRUE)
legend(125,-0.25,c("pdelta=0","0 < pdelta < 1","pdelta=1"),text.width=c(60,120,150),
col=c("#00ea6e","#D8D8D8","#ff0000"),lwd=2,lty=1,horiz=TRUE,bty="n")
## CI plot of difference
par(xpd=FALSE)
plot(time.diff,pred.diff,ylim=c(-0.1,0.35),xlab="",ylab="Difference",type="n")
axis(side=1)
mtext("Days",side=1,line=2.25)
rect(0,-0.05,max(time.diff),0.05,col=rgb(208,216,232,max=255),border=NA)
lines(time.diff,pred.2-pred.1,lwd=2,col="black")
abline(h=0,lty=2,lwd=0.5)
lines(time.diff,ub,lty=2,col="red")
lines(time.diff,lb,lty=2,col="red")
for (i in 1:length(time.diff)){rug(time.diff[i],col=COL[i],lwd=1.2,ticksize=0.04)}
axis(side=1)
par(xpd=TRUE)
legend(125,-0.21,c("pdelta=0","0 < pdelta < 1","pdelta=1"),text.width=c(60,120,150),
col=c("#00ea6e","#D8D8D8","#ff0000"),lwd=2,lty=1,horiz=TRUE,bty="n")
par(xpd=FALSE)
legend(675,0.37,c("Difference","95% CI"),col=c("black","red"),lty=c(1,2),lwd=2,bty="n")
###### Leukemia Example
data(leukstats)
hist(leukstats$t.stat, breaks=100, main="", xlab="")
#Plot sgpv no ordering
plotsgpv(est.lo=leukstats$ci.lo, est.hi=leukstats$ci.hi,
null.lo=-0.3, null.hi=0.3,
set.order=1:nrow(leukstats),
x.show=7000,
plot.axis=c("TRUE","FALSE"),
int.col="cornflowerblue",
null.pt=0, outline.zone=TRUE,
title.lab="Leukemia Example", y.lab="Fold Change (base 10)",
legend.on=FALSE)
axis(side=2,at=round(log(c(1/1000,1/100,1/10,1/2,1,2,10,100,1000),
base=10),2),labels=c("1/1000","1/100","1/10","1/2",1,2,10,100,1000),
las=2)
#Plot sgpv with ordering
plotsgpv(est.lo=leukstats$ci.lo, est.hi=leukstats$ci.hi,
null.lo=-0.3, null.hi=0.3,
set.order=order(leukstats$p.value),
x.show=7000,
plot.axis=c("TRUE","FALSE"),
null.pt=0, outline.zone=TRUE,
title.lab="Leukemia Example", y.lab="Fold Change (base 10)",
x.lab="Classical p-value ranking",
legend.on=TRUE)
axis(side=2,at=round(log(c(1/1000,1/100,1/10,1/2,1,2,10,100,1000),
base=10),2),labels=c("1/1000","1/100","1/10","1/2",1,2,10,100,1000),
las=2)
#Plot sgpv with ordering
plotsgpv(est.lo=leukstats$ci.lo, est.hi=leukstats$ci.hi,
null.lo=-0.3, null.hi=0.3,
set.order=order(leukstats$p.value),
x.show=1000,
plot.axis=c("TRUE","FALSE"),
null.pt=0, outline.zone=TRUE,
title.lab="Leukemia Example", y.lab="Fold Change (base 10)",
x.lab="Classical p-value ranking",
legend.on=FALSE)
axis(side=2,at=round(log(c(1/1000,1/100,1/10,1/2,1,2,10,100,1000),
base=10),2),labels=c("1/1000","1/100","1/10","1/2",1,2,10,100,1000),
las=2)
##Bonferroni
res = sgpvalue(est.lo=leukstats$ci.lo, est.hi=leukstats$ci.hi,
null.lo=-0.3, null.hi=0.3)
plot(c(1:nrow(leukstats)),
sort(leukstats$p.value),
ylim=c(0,1),
type="l",
main="Bonferroni adjusted vs. Raw p-values",
ylab="Probability",
xlab="Number of Rejected Hypotheses (index)")
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort( FDRestimation::p.fdr(leukstats$p.value)$fdrs), col="blue")
#Table of results
tab.res = table(res$p.delta>0, leukstats$p.value<(0.05/nrow(leukstats)))
colnames(tab.res) = c("SGPV=0", "SGPV>0")
rownames(tab.res) = c("p_Bon>0.05", "p_Bon<0.05")
kable(tab.res)%>%
kable_styling(c("striped", "bordered"))
| SGPV=0 | SGPV>0 | |
|---|---|---|
| p_Bon>0.05 | 65 | 164 |
| p_Bon<0.05 | 6799 | 100 |
#Put SGPVs on plot
plot(c(1:nrow(leukstats)),
sort(leukstats$p.value),
ylim=c(0,1),
type="l",
main=" Inference by Tail area, delta = 0.3",
ylab="Probability",
xlab="Number of Rejected Hypotheses (index)")
points(rank(leukstats$p.value,
ties.method = "first"),(res$p.delta), col="green",
cex=0.2)
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(leukstats$p.value), col="black")
lines(c(1:nrow(leukstats)), sort( FDRestimation::p.fdr(leukstats$p.value)$fdrs), col="blue")
#Put SGPVs on plot
plot(c(1:nrow(leukstats)),
sort(leukstats$p.value),
ylim=c(0,1),
type="l",
main=" Inference by Tail area, delta = 0.1",
ylab="Probability",
xlab="Number of Rejected Hypotheses (index)")
points(rank(leukstats$p.value,
ties.method = "first"), sgpvalue(est.lo=leukstats$ci.lo, est.hi=leukstats$ci.hi,
null.lo=-0.1, null.hi=0.1)$p.delta, col="green",
cex=0.2)
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(leukstats$p.value), col="black")
lines(c(1:nrow(leukstats)), sort( FDRestimation::p.fdr(leukstats$p.value)$fdrs), col="blue")
#Put SGPVs on plot
plot(c(1:nrow(leukstats)),
sort(leukstats$p.value),
ylim=c(0,1),
type="l",
main=" Inference by Tail area, delta = 0.05",
ylab="Probability",
xlab="Number of Rejected Hypotheses (index)")
points(rank(leukstats$p.value,
ties.method = "first"),
sgpvalue(est.lo=leukstats$ci.lo, est.hi=leukstats$ci.hi,
null.lo=-0.05, null.hi=0.05)$p.delta,
col="green",
cex=0.2)
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(leukstats$p.value), col="black")
lines(c(1:nrow(leukstats)), sort( FDRestimation::p.fdr(leukstats$p.value)$fdrs), col="blue")
#Put SGPVs on plot
plot(c(1:nrow(leukstats)),
sort(leukstats$p.value),
ylim=c(0,1),
type="l",
main=" Inference by Tail area, delta = 1e-06",
ylab="Probability",
xlab="Number of Rejected Hypotheses (index)")
points(rank(leukstats$p.value,
ties.method = "first"), sgpvalue(est.lo=leukstats$ci.lo, est.hi=leukstats$ci.hi,
null.lo=-0.000001, null.hi=0.000001)$p.delta, col="green",
cex=0.2)
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(pmin(leukstats$p.value*nrow(leukstats),1)), col="red")
lines(c(1:nrow(leukstats)), sort(leukstats$p.value), col="black")
lines(c(1:nrow(leukstats)), sort( FDRestimation::p.fdr(leukstats$p.value)$fdrs), col="blue")
International Consortium for Prostate Cancer Genetics (Schaid and Chang 2055; ICPCG 2018)
3,894 subjects: 2,511 cases & 1,383 controls
247,000 single-nucleotide polymorphisms (SNPs) from Chromosome 6
Goal: Identify interesting SNPs potentially associated with prostate cancer
https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000733.v1.p1
Access is controlled so we will simulate data here.
Here we have 10,000 genes with expression from 0 to 4. Patients 0-50 do not cancer and 51-100 do have cancer.
T-tests are conducted for each gene and p-values and CIs are taken from them.
Null Interval: [0.5, 1.3]
n=10000
pros.sim.data = matrix(NA, nrow=100, ncol=n)
pros.pvalue = NA
pros.lo = NA
pros.hi = NA
for(i in 1:n){
pros.sim.data[c(1:50),i] = runif(50, 0, 4)
pros.sim.data[c(51:100),i] = runif(50, 0, 1.5)
t.res = t.test(pros.sim.data[1:50,i], pros.sim.data[51:100,i])
pros.pvalue[i] = t.res$p.value
pros.lo[i] = t.res$conf.int[1]
pros.hi[i] = t.res$conf.int[2]
}
# Gene number on the x-axis, delta-gap on the y-axis, using an interval null hypothesis of
plotman(est.lo=pros.lo,
est.hi=pros.hi,
null.lo=0.5, null.hi=1.3,
set.order=NA,
type="delta-gap",
ref.lines=NA,
int.pch=16, int.cex=0.4,
title.lab="Prostate Simulated Example",
y.lab="Delta-gap",
x.lab="Position",
legend.on=TRUE)
# Gene number on the x-axis, -log10(classical p-value) on the y-axis.
plotman(est.lo=pros.lo,
est.hi=pros.hi,
null.lo=0.5, null.hi=1.3,
set.order=NA,
type="p-value",
p.values=-log10(pros.pvalue),
ref.lines=-log10(0.05),
int.pch=16, int.cex=0.4,
title.lab="Prostate Simulated Example",
y.lab=expression("-log"[10]*"(p-value)"),
x.lab="Position",
legend.on=TRUE)
# Second-generation p-value (SGPV) on the x-axis, -log10(classical p-value) on the y-axis
plotman(est.lo=pros.lo,
est.hi=pros.hi,
null.lo=0.5, null.hi=1.3,
set.order="sgpv",
type="comparison",
p.values=-log10(pros.pvalue),
ref.lines=c(-log10(0.05)),
int.pch=16, int.cex=0.4,
title.lab="Prostate Simulated Example",
y.lab=expression("-log"[10]*"(p-value)"),
x.lab="Second-generation p-value ranking",
legend.on=TRUE)
plotsgpv(est.lo=pros.lo,
est.hi=pros.hi,
null.lo=0.5, null.hi=1.3,
set.order=order(pros.pvalue),
x.show=10000,
plot.axis=c("TRUE","FALSE"),
null.pt=0, outline.zone=TRUE,
title.lab="Simulated Example", y.lab="Fold Change (base 10)",
x.lab="Classical p-value ranking",
legend.on=TRUE)
axis(side=2,at=round(log(c(1/1000,1/100,1/10,1/2,1,2,10,100,1000),
base=10),2),labels=c("1/1000","1/100","1/10","1/2",1,2,10,100,1000),
las=2)